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Abstract 

Hyperaccretion flows with mass accretion rates far above the Eddington rate have an N-shaped equilib¬ 
rium curve on the E-M plane (with S and M being surface density and mass accretion rate, respectively). 

The accretion flow on the lower E branch of the N-shape is optically thick, advection-dominated accretion 
flow (ADAF) while that on the upper one is neutrino-dominated accretion flow (NDAF). The middle branch 
has a negative slope on the E-M plane, meaning that the flow on this branch is secularly unstable. To 
investigate how the instability affects the flow structure and what observable signatures are produced, we 
study the time evolution of the unstable hyperaccretion flow in response to variable mass injection rates by 
solving the height-averaged equations for viscous accretion flows. When a transition occurs from the lower 
branch to the upper branch (or from the upper branch to the lower branch), the surface density rapidly 
increases (decreases) around that transition region, which induces locally enhanced mass flow (referred 
to as non-steady mass flow) into (out of) that region. We confirm that the non-steady flow can create a 
kind of disturbance and that it propagates over the whole disk. However, the non-steady mass flow is not 
strong enough to induce coherent transition over the whole disk, unless the mass injection rate varies with 
time. When the injection rate continuously changes, the neutrino luminosity varies intermittently, thus 
producing step-function-like light curves, as the radiation efficiency discontinuously changes every time the 
local transition occurs. The effects of changing the N-shape and possible observational consequences on 
the gamma-ray burst emission are briefly discussed in relation to gamma-ray bursts. 
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1. Introduction 

Gamma-ray bursts (GRBs) are the brightest explosions 
in the universe. They release energy up to ^ 10®^ ergs, 
which is on the order of 1 O“^M 0 C^, in the form of gamma- 
rays within a few seconds. From the spectral features 
and light curves of their prompt and afterglow emission, 
they are considered to be produced in ultrarelativistic jets. 
Although the picture of the central engine has not yet 
been established, several constraints are required from the 
observations. It should be able to launch ultrarelativistic 
jets, meaning that relativistic objects should be involved. 
It should have two different characteristic timescales: the 
variability timescale of their prompt emission (^ 10 “^ — 
10“^ s, e.g., Beloborodov et al. 2000; Ackermann et al. 
2010; Troja et al. 2015) and the duration timescale of the 
prompt emission (^ 2 — 1000 s). The latter timescale is 
much longer than the dynamical timescale of relativistic 
objects, but there is a way to explain all the requirements. 
The most promising model of the GRB central engine is 
a hyperaccretion flow (Narayan et al. 1992; Narayan et 


al. 2001). In this model an ultrarelativistic jet would be 
launched from a massive accretion disk around a stellar- 
mass black hole whose mass accretion rate is at most ~ 
0.01 — IMq s“^. Such a system is expected to form after 
a cataclysmic event such as the gravitational collapse of 
a massive star or the merger of a neutron star - neutron 
star binary or a black hole - neutron star binary (Piran 
1999). 

Since such an accretion flow is extremely optically 
thick with respect to photons, it cannot cool via elec¬ 
tromagnetic radiation. Instead, since the density and 
temperature of this accretion flow are very high (p ^ 
10^ g cm“^, T ^ 10^° K), it would cool via neutrinos, 
and such an accretion flow is often called as “neutrino- 
dominated accretion flows” (NDAF; Popham et al. 1999; 
Narayan et al. 2001; Di Matteo et al. 2002; Kohri & 
Mineshige 2002; Kohri et al. 2005; Gu et al. 2006; Ghen & 
Beloborodov 2007; Kawanaka & Mineshige 2007; Liu et al. 
2007; Kawanaka & Kohri 2012; Xue et al. 2013; see also 
section 10.6 of Kato et al. 2008). In this model, the vari¬ 
ability timescale can be interpreted as the dynamical (i.e.. 
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rotation) timescale of an accretion disk, while the burst 
duration timescale corresponds to the accretion timescale. 
As for how a hyperaccretion flow can produce an ultrarel- 
ativistic jet, two physical processes are mainly discussed 
in the literature: neutrino pair annihilation {vD —^ e“e^; 
Eichler et al. 1989; Asano & Fukuyama 2000; Zalamea 
& Beloborodov 2011; Suwa 2013) and magnetohydrody¬ 
namic (MHD) mechanisms such as the Blandford-Znajek 
(BZ) process (Blandford & Znajek 1977; McKinney & 
Gammie 2004; Hawley & Krolik 2006). 

In either scenario, it is implied that the energy deposi¬ 
tion rate (corresponding to the jet luminosity) as a func¬ 
tion of mass accretion rate has a step-function-like behav¬ 
ior at a certain critical mass accretion rate. According 
to the neutrino pair annihilation scenario, in fact, we ex¬ 
pect a big jump in the neutrino emission efficiency across 
the critical accretion rate, Mign ~ 0.003 — O.OIMq s“^ for 
fiducial parameters (see Kawanaka et al. 2013a). That 
is, when the mass accretion rate is above Mign tbe accre¬ 
tion flow becomes an efficient neutrino emitter (“NDAF” 
regime) and the energy deposition via neutrino annihila¬ 
tion is efficient, while when M is below Mign the accre¬ 
tion flow is in the “advection-dominated” (ADAF) regime 
and neutrino pair annihilation emission is not at all effi¬ 
cient. Similar behavior is expected in the MHD scenario, 
as well. Kawanaka et al. (2013a) investigated the MHD jet 
luminosity predicted from a hyperaccretion flow and found 
that it shows step-function-like behavior at M ^ Mign. In 
either case, therefore, one can expect a drastic change of 
jet luminosity from a hyperaccretion flow when its mass 
accretion rate attains the critical value ^ Mign • 

If there is any physical mechanism which can mod¬ 
ulate the mass accretion rate around Mign, the jet lu¬ 
minosity would change violently, which may explain to 
the highly variable GRB emission. Theoretically, sev¬ 
eral types of disk instabilities are known (Lightman 
& Eardley 1974; Shibazaki & Hdshi 1975; Shakura & 
Sunyaev 1976; Piran 1978), and the large-amplitude lu¬ 
minosity variations of accreting systems, such as those 
observed in dwarf novae and X-ray binaries are often 
ascribed to the accretion disk instabilities (e.g., Osaki 
1974; Mineshige & Wheeler 1989; Honma et al. 1991, see 
a concise review by Kato et al. 2008 chapter 5 and 10). 
However, little was known about the disk instabilities of 
a hyperaccretion flow, especially that with mass accre¬ 
tion rate around Mign. Recently, Kawanaka et al. (2013b) 
found that the thermal equilibrium curve at the innermost 
region of a hyperaccretion flow on the E-M plane has the 
structure like a character “N”, which has a negative slope 
at M ^ Mign- This means the existence of a secular (or 
viscous) instability of the accretion flow. 

In our study, we for the hrst time solve the time evo¬ 
lution of the unstable hyperaccretion flow and show that 
the luminosity evolution occurs intermittently as a conse¬ 
quence of a local transition which takes place in a narrow 
region. The plan of this paper is as follows: in Section 
2, we describe the equations for disk evolution and the 
modeled N-shaped thermal equilibrium curve. In Section 
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3 we show the results of our calculations with an emphasis 
on the consequence of the instability. Section 4 is devoted 
to the discussion on the dependence on the N-curve, the 
observable behavior, the comparison with the case of S- 
shaped curve, and future issues. 


2. Methods of Calculation 

2.1. Basic Equations for Viscous Diffusion 


In the present study, we solve the height-averaged basic 
equations for viscous diffusion of an axisymmetric accre¬ 
tion disk (see, e.g., Kato et al. 2008 chapter 3): 


dE _ 1 dM 
dt 2'Kr dr ’ 





(1) 

( 2 ) 


Here, E(r) = f p(r,z)dz (with p being mass density) is the 
surface density, r is the distance from the center of a black 
hole, M = —2'nrYiVr (with Vr being the radial velocity), v 
is the kinematic viscosity, and D is the angular velocity. 
Equation (1) describes mass conservation, while equation 
(2) represents angular momentum conservation. In our 
calculations, we adopt the pseudo-Newtonian force for a 
Kerr black hole (Artemova et al. 1996) for calculating D 
in equation (2); that is, from 


F = - 


We find 


GMbh 


r‘^-P(r — th)^ 


D = 


GMbh 


7-3 P 


( 3 ) 

( 4 ) 


where Mbh is the black hole mass, th is the radius of the 
event horizon, P = (rms/^’n) ~ 1 (with rms being the radius 
of the marginally stable orbit) is a constant to be specified 
later, and we assumed the balance between the gravita¬ 
tional force and the centrifugal force. Inserting equation 
(4) to equation (2), we obtain 


M = 67r- 


/ 

(r-rn) = 


2 (j 


[r - {I + l3)rYi]r 2 


dr 


/3 + 1 , Q_ a , 

r 2 (r-g^rn) 


/ -sttl 

(r-rn) 2 


i/E 


.(5) 


We next transform equations (1) and (5) to non- 
dimensional forms. For this purpose we introduce the 
following dimensionless variables: 


X = \/r/rout, (6) 

a^H = (7) 

(t = E/Eo, (8) 

T = t/to with = (9) 

m = M/MQ, ( 10 ) 

P = vYi/Mq. ( 11 ) 


Here, Tout is the size of the innermost region of the disk, 
which we are concerned with. Eg and Mg represent typical 
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surface density and mass accretion rate of NDAF, and to 
corresponds to the viscous timescale at the outer edge of 
the disk. These normalization constants will be specified 
later (see section 2.2). We can then rewrite the basic 
equations with dimensionless variables as follows: 

da 1 dm ^^2) 


dr A'kx^ dx 


{x^ — x^Y^^ d 
— (1 + d)xY\ dx 


m = 3TT- 


„/J+l - , 


l) 


{x^-xl 


13 + 2 


2. The entire equilibrium curve up-shifts as x increases. 
In other words, the bigger x is, the bigger is the value 
of fj, for a fixed a. 

3. The N-shape is more enhanced at smaller radii, 
whereas the unstable branch (with a negative slope) 
disappears at x > 0.5. 

In the present study, we set rout = 64i?g. Here, Rg{= 
GMbu/c^) is the gravitational radius. By comparing 


.(13)Figure 1(a) with Figure 1(b), 
ues to Mo and Eg: 


we assign the following val- 


2.2. The N-shaped Equilibrium Curve 

To solve the basic equations (12) and (13), we need to 
prescribe the functional form oi p = fJ.{a). This can be 
done by using the thermal equilibrium curve, on which 
heating and cooling are balanced. The thermal equilib¬ 
rium curve of the hyperaccretion flow has a characteristic 
N-shape, as is shown in the left panel of Figure 1 taken 
from Kawanaka et al. (2013b). They adopted a viscosity; 

i2=‘^aCsH (14) 

with a = 0.1, where H is the half thickness of the disk and 
Cs is the speed of sound. The vertical axis is M (mass ac¬ 
cretion rate) and the horizontal axis is E (surface density) 
in the left panel of Figure 1. We model this equilibrium 
curve in a simplified form as is shown in the right panel 
of Figure 1: 

For a < aA 


fj. = xa 


(15) 


For aA < a < as 

(x<0.5) 
(x>0.5) 

For a > aB 
a 

h- = Mb — 
o-B 


M = 


MA 

MA 




(16) 


(17) 


where a a and ctb are the values of cr at the two critical 
points, ‘A’ and ‘B’, respectively, where the slope of the 
thermal equilibrium curve changes; 


aA = 5.0-1- 


10.0 

075 


{x 


0.25), 


and aB = 15.0 


(18) 


Likewise, pLA and pB are the values of p. at the critical 
points, ‘A’ and ‘B’, respectively, and their explicit forms 
are calculated by using equations (15) and (16) as func¬ 
tions of X. 

Equations (15)-(17) represent the lower, middle, and 
upper branches of the thermal equilibrium curve, respec¬ 
tively. The principle for constructing this formula is as 
follows: 


Mo = 3.5 X 10-^[Mq s"^], (19) 

Eg = 5.5 X 10^^[g cm-2], (20) 

for the black hole mass of Mbh = 3Mq. Inserting these 
normalization constants into equation (9) and simply as¬ 
suming that to oc Mbh , we find 

to(MBH)=6.4xlO-2(^^)[s]. (21) 


2.3. Simplified Energy Equation 


As long as the system is in the thermal equilibrium, it 
evolves along the equilibrium curve. It starts to deviate 
from the curve, however, once the instability takes place. 
In order to describe the disk evolution in the state out 
of equilibrium, we introduce the following simple form of 
energy equation: 


dp, _ p- Peg 
dr Tth 


( 22 ) 


where peq is the value of p on the thermal equilib¬ 
rium curve for a given cr, and Tth represents the thermal 
timescale; 


rth — / * 'Tvis 


with 



, cr 4 

to-x 

M 


(23) 


Here, TvIs represents the diffusion timescale at radius r 
and is shorter at smaller r, and / is a constant set to be 
/ = 0.1 in the present study. 

2 . 4 . Boundary Conditions and Initial Conditions 


There are two boundary conditions. One is the choice 
of rhinj, mass injection rate at the outer edge of the disk 
(the outer boundary condition) and the other is the zero 
torque condition at the inner edge of the disk (the inner 
boundary condition); 

rh = mi„j(r) at r = Tout (24) 

p = 0 at r = Tin (25) 


In this study, we examine the following three cases for 
the way of mass injection from the outer boundary: 

Model 1 (slow rise): minj=TOgexp(r/10) withmg = 37rx8.60(26) 


1. The value of the kinematic viscosity is constant (i.e., 
p (X a) on the lower and upper branches, while p oc 
l/tT on the middle branch in the innermost region. 
The bigger the radius is, the smaller is the slope of 
the middle branch. 


Model 2 (rise): ?7iinj =TOgexp(r) with mg = 37rx8.60(27) 

Model 3 (decay): rhinj =rhgexp(— t) with mg = 37rx33.00(28) 

Note that the value of mg of Models 1 and 2 (or Model 
3) has been chosen in order that the whole disk can be 
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Fig. 1. Two thermal equilibrium curves: The left panel (a) the curve calculated by Kawanaka et al. (2013b). The right panel (b) 
the simplified curve used in the present study. The symbols ‘A’ and ‘B’ denote the two critical points. 


on the lower (upper) branch; i.e., rhinj < Stt x 8.69 (mi„j > 
Stt X 32.65). Also note that the rise timescale of Model 1 
is 10 to = 0.64 (Mbh/3Mq) [s] [see equation (21)]. 

The initial condition is a steady disk corresponding to 
toq. By integrating equation (13) with respect to a;, we 
obtain 




m 


3+2 




Stt 




3 - 1 : 

3 


)^H 


r/3 + 1 


r/5+1 




(x2-4)f (xl-xi)^ 


and the initial value of ^ is calculated by inserting mg into 
the equation above. We can then solve time-dependent 
equations (15)-(17) to obtain surface density evolution 
at each x. Actual numerical procedure is explained in 
Appendix 1. 

We set the inner boundary at Hn = 3.84i?g or Hn = 
0.24. The inner radius corresponds to the innermost stable 
circular orbit of a black hole with the dimensionless spin 
parameter a* of 0.6, for which th = 1.8i?g. We finally get 
P = (rms/rn) -1^ 1 . 12 . 

3. Results 

3.1. Basic Consideration Based on a Steady Model 

How does a transition from the lower branch to the up¬ 
per branch proceed in the disk when the mass injection 
rate gradually increases? To answer to this question, we 
first show in Figures 2(a) and 3(a) a series of the surface- 
density distribution expected to be the steady disk model 
for several mass injection rates; rhinj = rho exp(T/10) 
with T = 0.2, 0.8, 1.4, and 2.0 in Figure 2(a) while 
T = 0.0, 0.3, •••, 3.0 in Figure 3(a), respectively [recall 
equation (26)]. We, here, assume that the disk is on the 
lower (cr) branch of the N-shaped curve as long as cr < cta 
but that it is on the upper branch otherwise. Note that 
each curve in Figure 3(a) except for the lowest one is ver¬ 
tically offset by the value depending on the elapsed time. 

In these plots the mass accretion rates are taken to be 
the same everywhere by the assumption that every portion 
of the disk can immediately adjust to any changes in rhinj 
given at the outer boundary. As rhinj increases, rh at 


each radius gradually increases, so does the value of a. 
Eventually a transition from the lower ADAF branch to 
the upper NDAF branch takes place when cr > cta- We 
can see in Figures 2(a) and 3(a) that the surface density 
rapidly increases at the radius where an upward transition 
occurs. After the transition, neutrino cooling is efficient. 
Hence, neutrino is emitted from the high-cr regions. 
,(25|ie steady model is simple and easy to calculate, since 
it requires no numerical technique, however, there exist 
certainly limitations. For example, it should take time for 
the influence of varying mass injection rates to prevail over 
the disk plane, but such a time delay is totally neglected 
in the steady model. Also an upward transition cannot oc¬ 
cur in an infinitesimal time but a transition is assumed to 
be instantaneous in the steady model. Further, mass in¬ 
flow into and/or outflow from the transition region, which 
should be necessary to make a local change in cr, are not 
properly taken into account. We thus need to perform 
numerical simulations. 

3.2. Model 1 (slow rise) 

3.2.1. Surfaee-Density Evolution 

The first case is Model 1, in which the mass injec¬ 
tion rate slowly increases. Figures 2(b) and 3(b) display 
the global evolution of the surface-density distribution in 
slightly different ways. The expansion of the high-tr re¬ 
gion, where neutrino emission is significant, is clear in 
these figures. The upward transition to the upper (NDAF) 
branch is indicated by a rapid increase of cr. 

Let us carefully examine to what extent the steady 
model depicted in Figure 2(a) reproduces the time- 
dependent calculations shown in Figure 2(b). In the 
steady model, the mass accretion rate (rh) is everywhere 
equal to the mass injection rate (rhinj); in other words, 
the whole disk is assumed to immediately respond to 
a change of rhinj at the outer boundary. Another big 
assumption in the steady model is that a transition 
from one branch to another at a certain radius does 
not affect its environment. From these two panels, 
we notice that the time of an upward transition at a 
certain radius in the time-dependent model is delayed. 
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K/64/?g) 



K/64/?g) 


Fig. 2. Variety of surface-density (a) distributions based on the steady model (a) and on the time-dependent model (b) for various 
mass injection rates of rhinj = moexp(r/10) with r = 0.2, 0.8, 1.4, and 2.0. The critical lines (cta and ctb) are also plotted by the 
dashed line and the dash-dot line, respectively. 




Fig. 3. Series of a distributions in the steady model (panel a) and in Model 1 calculations (panel b). The corresponding mass 
injection rates are rfiinj = moexp(T/10) with r = 0.0, 0.3, ■■■, and 3.0. Each line is offset vertically by the value of 50/3 X r for 
visible clarify. 


compared with that at the same radius in the steady 
model. This feature can also be confirmed by comparing 
Figures 3(a) and (b) and is obviously the effect of the 
finite accretion timescale. We also find that the cr profile 
shows fluctuating patterns. These fluctuating patterns 
are more clearly displayed in Figure 3(b). This is caused 
by rapid mass exchanges between the neighboring regions 
associated with branch transitions. Hereafter we call such 
fluctuations as ‘non-steady’ effects. Obviously they are 
induced by ‘non-steady’ mass accretion flow (i.e., locally 
enhanced mass flow). We will discuss their physical 
meaning in details in the next subsection. 

3.2.2. Mechanism of Transitions Induced by Non-steady 
Effects 

In order to understand the non-steady effects, more 
careful analysis is necessary regarding the behavior of each 
part of the disk. This can be done by using the cr-fi plane 
at fixed radii. A typical example is displayed in Figure 
4. We see that the evolutionary path shown by the solid 
line grossly deviates from the thermal equilibrium curve 


4 


10 


0 





thermal equilibrium curve- 

time evolution - 


5 


10 ^ 


30 


log(a) 


Fig. 4. The actual behavior on the a-fi plane at the radius 
X = 0.39, where an upward transition is initiated at r < 0.357: 
The solid line represents the time evolutionary path, while 
the dashed line represents the thermal equilibrium curve. The 
arrow represents the direction of the time evolution. 


























































6 


Kimura et al. 


indicated by the dashed line. We can understand this fact 
in the following way. 

Suppose that the critical point A is reached at some 
radius. Although there are no longer low branch solutions 
above that point, a should steadily increase, since mass 
is continuously supplied from the neighboring outer zone. 
Thus, the disk should evolve in the right direction on the 
a-fj, plane. This is exactly what we see in Figure 4. 

After reaching the upper branch, however, unexpected 
behavior appears. That is, the evolutionary path fluctu¬ 
ates around the equilibrium point, on the upper branch. 
This is a result of complex mass exchanges between the 
neighboring sites. To understand how this occurs, we 
make the following thought experiments. 

At first we wish to note that the exact behavior of the 
disk on the a-fj, plane depends critically on the ratio of the 
thermal timescale (rth, on which /i changes) to the viscous 
timescale (tvIs, on which cr changes). Let us, hence, con¬ 
sider two hypothetical cases: Case A transition in which 
Tth/TVis <C 1 as is expected in a standard-type disk (see, 
e.g., Kato et al. 2008, chapter 3) and Case B transition in 
which Tth/rvis ~ 1. (The actual simulation will be shown 
later.) 

In Case A transition (see the upper four panels of Figure 
5), we assume rth/rvis "Cl. A short thermal timescale 
means that a system can quickly adjust to any changes 
added externally so that it should always stay around the 
thermal equilibrium state. [We assume that the disk is 
thermally stable (see, e.g., Kawanaka & Mineshige 2007), 
as in the case that we are discussing.] Let us denote by 
Tign the ignition radius, where the instability sets out. 
Since the disk should evolve along the thermal equilib¬ 
rium curve, <7 should increase accordingly, while /i should 
decrease during an upward transition at rjgn (see panel 
a2 of Figure 5). A non-steady mass flow (Am) is then 
induced from its neighboring sites (with high /r values) 
to the ignition region (with low p, values) as is illustrated 
in panel a4, since crudely we have rh oc d{x^)/dx from 
equation (13). Here, the non-steady mass flow is defined 
by 

Am(r) = m—(m), (30) 

where the averaged value, (rh), is the mass flow rate aver¬ 
aged over a region around the radius, r. This non-steady 
mass flow assures an increase of a at the ignition point, 
whereas a in the neighboring sites should decrease (see 
panels al and a3 of Figure 5). To conclude, a transition 
at a certain point does not trigger further transitions in 
the surrounding region. 

In order for a transition to propagate into the neigh¬ 
boring regions, these regions should acquire mass from 
the ignition site (at ngn)- For this to occur, the value of 

has to increase at rign, as is illustrated in panel b2 of 
Figure 5 (referred to as Case B transition); i.e., the disk 
evolutionary path on the plane should largely deviate 
from the equilibrium curve. For this sort of behavior to 
take place we should require Tth/rvis 1, but this never 
happens in a standard type disk with no instability, since 
roughly we estimate 


[Vol. , 

Tth/Tvis ^ (H/r)'^ < 1, (31) 

(with H being the scale-height of the disk), and since 
the aspect ratio is expected to be small, (H/r) <C 1. We 
wish to emphasize, however, that this inequality no longer 
holds, once an instability grows and produces abrupt, spa¬ 
tial variations of physical quantities. Under such circum¬ 
stances, equation (31) should be written as 

Tth/Tvis (H/Ar)^ ^ 1, (32) 

with Ar being the spatial scale, on which the value of /i 
varies. If this were the case, the evolution of a profile 
would be like that illustrated in panel b4 of Figure 5. 

Which is actually the case in the time-dependent calcu¬ 
lation? The resultant time evolution of an unstable disk is 
illustrated as Case C in the lower four panels of Figure 5. 
When a transition occurs, /i only slightly decreases at ngn 
[see Figure 4]. This apparently looks like Case A except 
for one important difference; that is, the disk trajectory 
starts to oscillate around its equilibrium value through in¬ 
teractions with neighboring regions after the transition to 
the upper branch. (The meaning of ‘interactions’ will be 
discussed in more details in Sec. 3.2.3.) Owing to these 
fluctuations, the values of ^ in the surrounding regions 
also begin to oscillate around their equilibrium values, and 
such fluctuations in ^ eventually induce mass accretion 
flow from or into the neighboring region, thereby occa¬ 
sionally triggering transitions [see Figure 5(c)]. These are 
purely non-steady effects. 

Importantly, the onset of the instability at some radius 
does not always trigger the instability in adjacent regions. 
The instability zone is conhned in a relatively narrow 
region. In other words, the instability does not promote 
coherent behavior of the entire disk. This point will be 
discussed in section 4.4 in comparison with the case of 
disks with the ‘S’-shaped equilibrium curve. 

3.2.3. Wave Propagation 

There is another interesting feature unique to unsta¬ 
ble accretion disks; the influence of rapid variations in rh 
propagates from rign in the radial direction, as is depicted 
in Figure 6(a). We can see a wavy pattern in the rh dis¬ 
tribution propagating over the disk plane from rign both 
in the inner and outer directions. Fluctuations in m drive 
small amplitude fluctuations in a [see equation (12)]. 

Figure 6(b) illustrates the time developments of the ac¬ 
cretion rate (to) at several fixed radii. Here, by m at r 
(= ngnTigniAr, ■■■) we mean the mass flow from the 
radius r 4- Ar to the radius r. The upper three curves (or 
the lower three curves) represent non-steady mass inflow 
(outflow) towards the ignition point from the outer (in¬ 
ner) adjacent zones. It is important to note time delays 
in the to response at the distant points (in comparison 
with that at the closer points) to the ignition radius, in¬ 
dicating that the fluctuations are not locally conhned but 
spatially propagate and cause cr variations in the neigh¬ 
boring sites. When the values of cr exceed its critical value, 
(ta, an upward transition can be triggered. The transition 
can occur, even when a is expected to be below cta in the 
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Fig. 5. Schematic pictures explaining how a transition occurs by non-steady effects. Here, we consider two hypothetical cases in 
the upper and middle panels and one actual calculation results in the bottom panels. Left nine panels illustrate time-dependent 
evolution of the disk on the a-fi plane and associated mass flow at the ignition point rjgn (panels a2, b2, and c2) and at its neighboring 
sites (other panels). Note that Am (non-steady mass flow) represents mass flow other than the global mass flow (see text for the 
definition). Right three panels a4, b4, and c4 illustrate the time-dependent evolution of the disk on the r-a plane. 




Fig. 6. The vertical axis “mass accretion rate” represents rh. (a) Propagation of the m variation wave due to the instability over 
the disk plane. The elapsed times are 0.334, 0.337, 0.340, •••, and 0.355, and for visible clarity each curve is offset upward by 0, 12, 
24, •••, and 84, respectively. The two arrows indicate the times and places of the onset of the instability, (b) Time variations of rh 
at several radii illustrated on the r-m plane. Each curve is offset by —30, — 20, — 10, • • •, and 20, respectively. 
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steady model for a given minj. 

Here, we discuss what is meant by “interactions” in 
Sec. 3.2.2. The unique features in Case C transitions 
reside in back-reactions of the non-steady mass flow 
(Am) [see the upper and lower panels of Figure 5]. 
According to our simulations, a and /r values change in 
the neighboring region after their receiving non-steady 
mass inflow. The changes in ^ will then modify the mass 
flow pattern in a next moment, since m oc d{x^)/dx. The 
modified mass flow pattern in the next moment produces 
further changes in a and /i values, which sometimes cause 
a negative feedback. This is the back-reaction of the 
non-steady mass flow, which is not considered in Case A 
transition. The transition wave eventually damps but on 
longer timescale than the thermal time. This explains 
why the solution oscillates around the equilibrium curve 
even after the thermal timescale. 


3.2.4- Neutrino Luminosity and Mass Accretion Rate 
The neutrino luminosity can, in principle, be calculated 

by 

/•i 


with 
Lo = 


r™* GMbhM 
Irin 2r^“^(r — th)^ 

GMbhMq 

Rz 


dr = L 


m 


° 64a;^ 

(34) 


Cautions should be taken here, since m in equation (33) 
should be mass flow rate on the upper (NDAF) branch. 
We, hence, need special arrangements when calculating 
neutrino luminosity from the part on the lower branch, 
which has very little contribution to the neutrino lumi¬ 
nosity, even if m is relatively large. In actual calculations, 
therefore, we use equation (17) even when ct < ctb to cal¬ 
culate /r from cr, and then calculate m from equation (29) 
by assuming the steady-disk relation. 

The resultant neutrino light curve is shown in Figure 
7(a). In the same figure, we plot two more light curves rep¬ 
resenting the cases with no N-shaped equilibrium curves 
nor instability. The upper (or lower) straight line la¬ 
beled with ‘steady l’(‘steady 2’) is the light curve that 
we would expect if the lower (upper) branch would ex¬ 
tend up (down) to high (low) a regimes. 

Comparing these light curves, we find two key features. 
First, the overall shape of the light curve is steeper than 
the cases with no instability. The main reason for this 
resides in the ‘N’-shape; that is, the value of a and, hence, 
the /i-value suddenly increases on a shorter timescale than 
the viscous timescale (on which rhinj increases) after the 
onset of the instability. 

The second key feature of the light curve is that the lu¬ 
minosity abruptly increases every time a transition takes 
place at a certain radius, resulting in step-function-like 
variations. This is because that the instability is confined 
in a narrow region and does not produce coherent varia¬ 
tions over a wide range of the disk plane, as mentioned 
above. 
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We compare with the mass accretion rate at radius Tin 
with the mass injection rate in Figure 7(b). We can see 
significant variability in the mass accretion rate. This may 
produce rapid variations in the electromagnetic radiation 
via the Blandford-Znajek process (to be discussed in sec¬ 
tion 4.3). 

3.3. Model 2 (rise) 

The second case is Model 2, in which the mass injection 
rate increases more rapidly. We show the global evolution 
of the surface-density distribution in Figure 8(a). We can 
see that the region on the upper (NDAF) branch expands 
from rign both in the inner and the outer directions more 
rapidly than in Model 1. This is simply because of more 
rapid increase in rhinj in Model 2 than in Model I. If we 
have a closer look at this figure, moreover, we see that 
the out-going transition wave propagates faster than the 
in-going wave. 

We show in Figure 8(b) the neutrino luminosity in 
Model 2. We can see that the luminosity rises more 
smoothly in Model 2 than in Model I because the speed 
of the instability propagation is faster than in Model I. 

^^^^^hodel 3 (decay) 

3 . 4 . 1 . Surface-Density Evolution 

In the last model (Model 3), we let the mass injection 
rate decrease on the viscous timescale. We compare 
the steady model (left) and the time-dependent model 
(right) in the same way as we did in Model I [see Figure 
9]. The adopted mass injection rates are mo exp(—r) 
with T = 0.0, 0.8, 1.4, and 2.0 [see equation (28)]. We 
clearly see in both these panels that the high-cr region 
producing intense neutrino emission gradually shrinks. 
In addition, we notice a tendency that the high-cr regions 
are systematically narrower in the steady model than in 
the time-dependent model. This is because it takes about 
the viscous timescale until the information of reduced 
rhinj is transported to each portion of the disk so that it 
can lose the surface density to reach point B in the right 
panel of Figure I. 

3 . 4 . 2 . Mechanism of Transitions Induced by Non-steady 
Effects 

Figure 10 shows the simulation results of Model 3 on 
the a-yi plane, while Figure II is a schematic picture 
explaining the time evolution of each portion of the disk 
on same plane. When point B is reached at some radius 
(referred to as Xign), a downward transition is triggered 
there, thereby the value of yi sightly increasing. This 
then causes non-steady mass flow (Am) from rjgn to 
the neighboring regions so that cr should decrease at 
Tign. Even after the downward transition is completed, 
the /c and cr are oscillating around their equilibrium 
values more remarkably than in Model I. We explain 
why these violent oscillations occur in Sec. 3.4.3. Such 
fluctuations in /i induce mass accretion flow from or 
into the neighboring regions and occasionally trigger 
transitions to the lower branch in the surrounding region. 
It is also important to note that such propagation of the 





Time Evolution of a unstable NDAF 


9 


No. ] 




Fig. 7. (a) The solid line represents the calculated neutrino luminosity in Model 1. The two dashed lines (‘steady 1’ and ‘steady 2’) 
show what the light curves would be in case that the thermal equilibrium curve were not N-shaped; that is, we assume in ‘steady 1’ 
that the lower branch [equation (15)] would extend up to the higher a regime while in ‘steady 2’ the upper branch [equation (17)] 
would extend down to the lower a regime, (b) Time variations of mass accretion rate at the inner edge of the disk (solid line) and 
of the mass injection rate (dashed line), respectively. 




Fig. 8. (a) Same as Figure 2(b) but for Model 2. The corresponding mass injection rates are minj =moexp(r) with r = 0.2, 0.8, 1.4, 
and 2.0. (b) Same as Figure 7(a) but for Model 2. 




Fig. 9. Same as Figure 2 but for Model 3 (decay). We compare the results of (a) the steady model and (b) the time-dependent 
model at the elapsed times of r = 0.0, 0.8, 1.4, and 2.0. 
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Fig. 10. The actual behavior on the a-fi plane at the radius 
x = 0.27 in Model 3: The solid line represents the time evo¬ 
lutionary path, while the dashed line represents the thermal 
equilibrium curve. The arrow indicates the direction of time 
evolution. 

instability is restricted in a narrow zone. In other words, 
it cannot produce coherent oscillation. 

3.4- 3. Wave Propagation 

Next, we show in Figure 12 the radial distribution of 
the mass accretion rate after the downward transition. It 
is clear in this figure that the rh distribution is never flat, 
meaning that each part of the disk at any times suffers 
fluctuations. Especially, m rapidly goes up and down with 
a large amplitude when a transition takes place. The in¬ 
fluence propagates in a form of m variation wave from rign 
to both in the inner and outer directions. The situations 
are similar to those in Model I but are more pronounced 
in Model 3. Fluctuations in m drive variations of a [see 
equation (12)]. 

Here, we think the violent fluctuations mentioned in 
Sec. 3.4.2 over again [see Figure II]. The fluctuations 
are the results of the back-reactions of the non-steady 
mass flow as mentioned in Model I. The values of cr and 
p, change in the neighboring zones after their receiving 
non-steady mass outflow. The changes in p will then 
modify the mass flow pattern in a next moment. The 
modified mass flow pattern in the next moment further 
changes cr and p values, which sometimes cause a negative 
feedback. For this reason, the a and p values violently 
oscillates around the equilibrium curve even after the 
thermal timescale. 

3.4.4- Neutrino Luminosity and Mass Accretion Rate 
We show in Figure 13(a) the neutrino luminosity in 

Model 3. Comparing the changes of the time-dependent 
luminosity with the changes of the luminosities in the 
steady models, we can see that the luminosity rapidly de¬ 
creases when downward transitions from the upper regime 
(NDAF) to the lower (ADAF) branch take place. 

We plot the accretion rate at the inner edge of the disk 
for Model 3 in Figure 13(b). We again see the significant 
variability in the mass accretion rate, as we saw in Model 


Fig. 12. Same as Figure 6(a) but for Model 3. The elapsed 
times are 1.711, 1.714, 1.717, and 1.732. Note that the 
curve is offset by 0, 12, 24, • • •, 84 for visible clarity. 

1 . 

4. Brief Summary and Discussion 

4-1. Brief Summary 

The thermal structure of the hyperaccretion flows is 
characterized by an N-shaped thermal equilibrium curve: 
the middle branch, which is secularly unstable, is sand¬ 
wiched by the lower ADAF branch and the upper NDAF 
branch. In this study, we examine the time evolution of a 
disk, which has N-shaped thermal equilibrium curves, by 
means of simple numerical simulations. We solve time- 
dependent, height-averaged basic equations for viscous 
disks receiving with variable mass injection to see how the 
transition between the lower and upper branches occurs, 
how it propagates over a disk plane, and what differences 
arise from the cases without N-shaped equilibrium curves. 
We studied three models, which have different time de¬ 
pendence of the mass injection rate: Model 1 (slow rise). 
Model 2 (rise), and Model 3 (decay). What we find can 
be summarized as follows: 

1. When the mass injection rate increases (or de¬ 
creases) with time, the region initially on the 
lower (upper) branch successively undergoes upward 
(downward) transitions. The transition associated 
with rapid changes in the surface density and prop¬ 
agates over the disk, giving rise to an expansion 
(shrinkage) of the neutrino-emission region. 

2. Non-steady mass flow (Am), the mass flow deviated 
from the spatially averaged flow, arises in a narrow 
region around the ignition radius rign during a tran¬ 
sition. This is because the middle branch is secularly 
unstable; that is, the larger (or the smaller) the sur¬ 
face density of the unstable region is, the more (less) 
mass enters the ignition region than its neighboring 
regions. Resultant m variation wave propagates over 
the whole disk. 

3. Such non-steady effects can occasionally cause tran¬ 
sitions even in the region at which the transition 
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Fig. 11. Same as Figure 5(c) but for Model 3 (decay). 



Fig. 13. Same as Figure 7{: 

would not occur, if the disk were in a steady state. 
The more slowly the mass injection rate changes, 
the more remarkable the non-steady effects become. 
The effects are, however, not strong enough to pro¬ 
duce coherent transition over the entire disk plane. 

4. The transition between the lower and upper 
branches yields abrupt changes in the neutrino lu¬ 
minosity. In fact, the neutrino luminosity calculated 
based on the time-dependent results changes much 
more rapidly than that in the case without N-shaped 
equilibrium curves. Non-steady effects can be seen 
as fluctuations of the neutrino luminosity and in 
mass accretion rates. 

4-2. Effects of the modified N-shape on disk evolution 

We have argued that the instability is a source of pro¬ 
ducing the non-steady effects. If so, the properties of the 
disk evolution may critically depend on the precise shape 
of the N-shaped thermal equilibrium curve. In this sub¬ 
section, we modify the thermal equilibrium curve in such 
a way to enhance the instability to see if it is possible for 
the instability to propagate more widely on the disk plane. 

Kawanaka et al. (2013b) adopted the a viscosity pre¬ 
scription in calculating the thermal equilibrium curve of 
a hyper accretion shown in the left panel of Figure 1. 
Although they set a to be constant (= 0.1), there are 
no strong reasons to believe so. As a thought experiment 
we set a = 0.01 on the upper branch, keeping the location 



t 


i) and (b) but for Model 3. 



Fig. 14. Modification of the N-shaped equilibrium curve at 
r = 4i?g: The solid line is the original curve while the dashed 
line is the modified one. 

of the lower branch and the slope of the middle branch un- 
changed. The value of a-Q is now 45.0, instead of 15.0 (see 
Figure 14). By this change, the middle unstable branch 
in the equilibrium curve is more elongated than in the 
original model. This model is referred to as the modified 
N-shape model. 

As to the mass injection rate, we only calculate the case 
of slow rise [i.e.. Model 1 in section 3; see also equation 
(26)]. First, we show the global evolution of the surface- 
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Fig. 16. Same as Figure 7(a) and (b) 



Fig. 15. Same as Figure 3(b) but for the modified N-shape 
model. The elapsed times are r = 0.0 — 3.0 with a constant 
interval of 5t = 0.3. Each line is offset vertically by the value 
of 100 X r for visible clarify. 

density distribution in Figure 15. We can see that the re¬ 
gion expands in a similar way to that in Figure 3(b). Note 
also that the high-cr region with efficient neutrino emission 
values in the upper branch are higher. Next, we calculate 
the luminosity and show it in Figure 16(a). Comparing 
it with Figure 7(a), we understand that the modification 
of the N-shape does yield more rapid increase in neutrino 
luminosity, but the speed of the upward transition wave 
is not so much affected. We also calculate the mass ac¬ 
cretion rate at radius Hn in Figure 16(b). Comparing it 
with Figure 7(b), we find that the accretion rate varies 
more or less continuously. These variations look similar 
to the observed ones in gamma-ray bursts (see the next 
subsection). 

4-3. Implications on the GRB Emission Feature 

Here we briefly discuss the observable features of the 
unstable hyperaccretion flows. The short-term variability 
of a hyperaccretion flow would be the origin of the highly 
variable prompt emissions of GRBs because the variable 
energy deposition from an accretion flow would produce 
a spatially inhomogeneous outflow, and it would result in 
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but for the modified N-shape model. 

producing internal shocks, which give rise to the highly 
variable emission from shock-accelerated electrons. 

According to our calculations, the neutrino luminosity 
shows abrupt changes with time during the transition of 
a hyperaccretion flow. In Model 1 with slowly increasing 
mass injection, especially, the step-function-like changes 
in the neutrino luminosity occur repeatedly. If we adopt 
the neutrino pair annihilation scenario as the energy depo¬ 
sition process from a hyperaccretion flow, we can expect 
the intermittent energy ejection from an accretion flow, 
which can result in the highly variable emission of GRBs. 

As for Model 3, where the injection rate decreases with 
time, the mass accretion rate at the inner edge shows 
strong variability. This may be related to the short-term 
variation of the Blandford-Znajek jet luminosity since it is 
proportional to the magnetic energy density on the event 
horizon and is, hence, likely to be proportional to the mass 
accretion rate (Kawanaka et al. 2013a). 

We wish to note that the behavior of the secular insta¬ 
bility is not always controlled by the variation timescale 
of the mass injection rate, but that rapid variations can 
occur on a much shorter timescale, less than 0.1 to ~ 10“^ 
s [see figure 7(b) and equation (21)]. We here need to 
distinguish two timescales: the timescale of variable mass 
injection rate (which is roughly the viscous timescale at 
large radii) and that of the fluctuations in the cr and p 
values around the ignition point (which is much shorter 
than the former). The variability timescale observed in 
GRBs is the order of ^ 10“^ s (e.g., Ackermann et al. 
2010; Troja et al. 2015). The observed rapid variability 
of the GRBs should be related to the latter process, and 
so we conclude that the variability timescale in our model 
can match the observations of GRBs. 

4-4- Comparison with the Dwarf-Nova Type Instability 

It is well known that the outbursts of dwarf-novae, 
a subclass of cataclysmic variables, and probably X-ray 
nova eruptions are triggered by a disk instability caused 
by the partial ionization of hydrogen in the disk (see, 
e.g., Kato et al. 2008, chapter 5 and references therein). 
This is the so-called dwarf-nova type instability (or hydro¬ 
gen ionization instability) and is characterized by the S- 
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shaped equilibrium curve. It can be easily shown that the 
middle branch is thermally and secularly unstable, while 
the other two branches are stable. Under the circum¬ 
stances relaxation oscillations between the bi-modal stable 
branches spontaneously occur, even if the mass injection 
rate does not vary in time. This is the most successful 
time-dependent theory of the viscous accretion disk. 

There are similarities and differences between the 
dwarf-nova type instability and the present case of hyper¬ 
accretion disk instability with the N-shaped equilibrium 
curve. The most prominent common feature is bi-modal 
nature of the equilibrium curve. There are three states for 
a given minj and the middle state being secularly unstable 
in the present case, while dwarf-nova type disks undergo 
bi-modal transitions, although only one state is found for 
a given rhjnj. Another common feature is global propaga¬ 
tion of bi-modal transitions from the ignition point to its 
neighboring sites. 

A big distinction resides in the fact that the solution 
is only secularly unstable in the present case. This is a 
crucial point when discussing the global nature of the in¬ 
stability in the present case, since the thermal timescale is 
much shorter than the viscous timescale, on which a secu¬ 
lar instability grows. We in fact find non-steady features 
not so much pronounced, judging from the fact that the 
transition is restricted in a narrow region and light curves 
exhibit small-amplitude fluctuations. This contrasts the 
case in the dwarf-nova type disks which undergoes co¬ 
herent oscillation and produce large-amplitude luminosity 
variations. To summarize, the big distinction arises from 
the fact that the /i value can change with a large ampli¬ 
tude, thereby producing significant non-steady effects in 
a thermally unstable disk, while it cannot in a thermally 
stable disk. 

4-5. Study in the Future 

It might be important in this respect that we have em¬ 
ployed one-dimensional formulation for the basic equa¬ 
tions in the present study, assuming that the vertical 
structure of the disk can immediately respond to changes 
in (T. This is not precisely correct, since it takes about a 
sound crossing time over the disk height (~ iJ/cg) for the 
entire vertical structure to adjust to any changes in a. In 
future work, we thus need two-dimensional simulations to 
see vertical propagation of the instability. 

Further, discontinuous transitions in the neutrino light 
curve, which seems to depend on a finite mesh size, can be 
altered by considering the two-dimensional effects. This is 
because changes in the vertical structure take longer time 
than the propagation time of the instability in the radial 
direction over a distance of one mesh size. We naively 
expect that the region with the width oi ^ H may simul¬ 
taneously undergo the branch transition, if we properly 
consider the two-dimensional effects. 

In order to solve the time evolution of a hyperaccretion 
disk strictly, we have to solve equation (22) without mod¬ 
eling it. In other words, we have to evaluate advection 
cooling rate as the function of the radial velocity 

(vr) and the r differentiation of the entropy (ds/dr) and 


neutrino cooling rate {Q^ ) as the function of the temper¬ 
ature (T). 

We are grateful to the anonymous referee for insightful 
comments. This work is supported in part by Grants-in- 
Aid of the Ministry of Education, Culture, Sports, Science 
and Technology (MEXT) (26400229 SM). 


Appendix 1. Numerical Procedures 

In our calculations, we use the implicit method for time 
integration of equations (12) and (13). At each time step, 
we solve these equations, as well as the energy equation 
in non-dimensional forms [see equation (22)]. The radial 
coordinate {x= 0.0 - 1.0) is divided into 200 mesh points 
(with a constant mesh spacing of Ax = 0.005); Xi = {Ax)i 
with i = 1,2, • • • ,200. We distinguish cell variables (such as 
a and fi) and interface variables (such as m). The former 
variables are defined at the same place of the radial mesh 
points; e.g., ai is the surface density at Xi (i = 1,2, ...200), 
while the latter variables are defined between the mesh 
points; e.g., m is the mass flow rate from the cell at Xi^i 
to the cell at Xi. 

Second, we set At to be constant (= 10“®) being inde¬ 
pendent of X. In the present study we assign At = 10“® 
so as to satisfy the two conditions, Acrjcr < 0.03 and 
A/r//r < 0.03 at every radius. When either of these condi¬ 
tions is unsatisfied at some mesh point (s), we re-calculate 
by using smaller At so that conditions should be met. 
One exception is Model 3, in which we omit calculations 
in the innermost region (at x < 0.26) to avoid numerical 
difficulties. 

Appendix 2. Accuracy of Numerical Calculatiou 

As a test of the numerical code, we check the mass 
conservation, we compute M^isk, the disk mass in the fol¬ 
lowing two methods and the difference between them. We 
use trapezoidal rule for this integral calculus. In method 
A we calculate 

pUout P^out 

= 27r / radr = Att j x^adx. (Al) 

J Uin Xi-n 

In method B we calculate 

-^disk= / {mni-rh^cc)dT. (A2) 

We have confirmed that the total mass is conserved 
within the accuracy of less than 0.2% in Model 1, less 
than 0.8% in Model 2, and less than 3.6% in Model 3 for 
time integration of r = 0 — 2.0 (rf = 2.0). 
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